# Import original data
d = import(here("03_updated_analyses_who", "data", "oxygen_and_corticosteroids_who_ma.xlsx")) %>%
# Only mortality outcome and patients using corticosteroids
filter(outcome == "mortality",
corticoid == 1)
# Find what studies had 0 patients in at least one of the treatment arms per
# subgroup
index_low =
d %>%
filter(trt_total == 0 | control_total == 0) %>%
filter(oxygen == "low")
index_niv =
d %>%
filter(trt_total == 0 | control_total == 0) %>%
filter(oxygen == "NIV")
index_imv =
d %>%
filter(trt_total == 0 | control_total == 0) %>%
filter(oxygen == "IMV")
# Only include studies that do not have 0 patients in a treatment arm (per subgroup)
d_clean =
d %>%
filter(case_when(
oxygen == "low" ~ study %nin% index_low$study,
oxygen == "NIV" ~ study %nin% index_niv$study,
oxygen == "IMV" ~ study %nin% index_imv$study)
)
# Calculate log odds ratio
d_logOR =
escalc(
measure = "OR", # log odds ratio,
# Tocilizumab
ai = trt_events,
n1i = trt_total,
# Control
ci = control_events,
n2i = control_total,
data = d_clean
) %>%
as_tibble() %>%
# Calculate standard error
mutate(sei = sqrt(vi),
# Set order to use "IMV" as the Intercept in brms
oxygen = factor(oxygen,
levels = c("IMV",
"low",
"NIV")))
# saveRDS(d_logOR,
# here("03_updated_analyses_who", "output", "data", "effect_sizes.Rds"))
The Bayesian meta-regression random-effect model is:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu & = \alpha_{subgroup[k]}\\ \\ \alpha & \sim \operatorname{Normal}(0, 0.82^2) \tag{Priors} \\ \tau & \sim \operatorname{Half-Normal}(0.5^2) \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio of tocilizumab versus control and \(\sigma_i^2\) is the known sampling variance in study \(i\). Because this is a random-effect model, each study \(i\) has its own distribution, where \(\theta_i\) represents its mean effect. All \(\theta_i\)s are drawn from normal distribution where the mean effect is \(\mu\) and the variance \(\tau^2\), which represents the between-study heterogeneity. \(\mu\) is predicted by a no-intercept linear regression, where each subgroup \(k\) (simple oxygen only; noninvasive ventilation; invasive mechanical ventilation) has its own parameter \(\alpha\), representing the mean effect of each respective subgroup.
As suggested by Depaoli and van de Schoot, 2017:
Our goal is to use priors that exclude impossible treatment effects, and yet applies little influence in the final results (hereafter, weakly informative priors).
We find highly unlikely that a pharmacological treatment, such as tocilizumab, will yield a 80% odds reduction in 28-days all-cause mortality. Thus, regarding \(\alpha\), we set a prior which the 95% CI ranges from 0.2 to 5.0 odds ratio, i.e. \(Normal(0, 0.82)\) in the log odds ratio scale.
sd_alpha = 0.82
twosd_alpha = sd_alpha*1.96
# https://stackoverflow.com/questions/19950219/using-legend-with-stat-function-in-ggplot2
ggplot(data = data.frame(x = c(-3, 3)), aes(x)) + #Empty plot
# Normal(0, 0.82)
stat_function(fun = dnorm, n = 1000,
args = list(mean = 0, sd = sd_alpha), linetype=1, size = 1.2,
aes(colour = "0.82")) +
scale_colour_manual(latex2exp::TeX("Normal(0, $\\sigma)"),
values = pnw_palette("Starfish", 3, type = "discrete")) +
geom_vline(xintercept = c(-twosd_alpha, twosd_alpha),
linetype = 2,
color = pnw_palette("Starfish", 3, type = "discrete")[1]) +
scale_y_continuous(breaks = seq(0, 0.6, 0.3),
limits = c(0, 0.6),
expand = c(0, 0)) + # remove gap between X and Y axis
scale_x_continuous(breaks = c(-twosd_alpha, 0, twosd_alpha),
labels = function(x) round(as.numeric(x), 2),
expand = c(0, 0)) +
coord_cartesian(x = c(-3, 3)) +
labs(x = "\nLog Odds Ratio",
y = "Density\n") +
theme_classic() +
theme(
plot.margin = margin(20,20,0,20),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
legend.position = 'right',
legend.text = element_text(size=12),
legend.title = element_text(size=14),
legend.key= element_blank(),
panel.background = element_blank()
)
95% quantile intervals:
tibble(
Mean = 0,
"SD" = sd_alpha,
) %>%
mutate("Mean " = exp(Mean),
"95% CI" =
str_c(" [",
round(exp(Mean - 1.96*`SD`), 1),
", ", round(exp(Mean + 1.96*`SD`), 1),
"]")) %>%
gt() %>%
cols_align(
align = "center",
columns = everything()
) %>%
# Add column tabs based on the column names defined above
tab_spanner(label = "Log Odds Ratio",
columns = c("Mean", "SD")) %>%
tab_spanner(label = "Odds Ratio",
columns = c("Mean ", "95% CI"))
| Log Odds Ratio | Odds Ratio | ||
|---|---|---|---|
| Mean | SD | Mean | 95% CI |
| 0 | 0.82 | 1 | [0.2, 5] |
In the log odds ratio scale, Spiegelhalter et al. 2004 categorized \(0.1\) - \(0.5\) as “reasonable”, \(0.5\) - \(1.0\) as “fairly high”, and \(> 1.0\) as “fairly extreme” heterogeneity ranges. Below, we show the percentages of density of each prior distribution of \(\tau\) for Spiegelhalter’s ranges. We also added a “low” category for values between \(0\) and \(0.1\).
The \(\operatorname{Half-Normal}(0.5)\) distribution yields plausible probabilities in each of these ranges.
sd_w = 0.5
median_w = round(extraDistr::qhnorm(0.975, sigma = sd_w),2)
# https://stackoverflow.com/questions/19950219/using-legend-with-stat-function-in-ggplot2
ggplot(data = data.frame(x = c(0, 2)), aes(x)) + #Empty plot
# Half-Normal(0.5)
stat_function(fun = extraDistr::dhnorm, n = 1000,
args = list(sigma = sd_w), linetype=1, size = 1.2,
aes(colour = "0.5")) +
scale_colour_manual(latex2exp::TeX("Half-Normal($\\sigma)"),
values = pnw_palette("Bay", 1, type = "continuous")) +
scale_y_continuous(breaks = seq(0, 2, 1),
limits = c(0, 2),
expand = c(0, 0)) + # remove gap between X and Y axis
scale_x_continuous(breaks = c(0, 0.1,median_w, 0.5, 1, 2),
labels = function(x) round(as.numeric(x), 1),
expand = c(0, 0)) +
geom_vline(xintercept = median_w, linetype = 2) +
coord_cartesian(x = c(0, 2)) +
labs(x = "\nLog Odds Ratio",
y = "Density\n") +
theme_classic() +
theme(
plot.margin = margin(20,20,0,20),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
legend.position = 'right',
legend.text = element_text(size=12),
legend.title = element_text(size=14),
legend.key= element_blank(),
panel.background = element_blank()
)
tibble(
"Low" = c(
100*round(phnorm(0.1, sigma = sd_w) - phnorm(0, sigma = sd_w),2)
),
"Reasonable" = c(
100*round(phnorm(0.5, sigma = sd_w) - phnorm(0.1, sigma = sd_w),2)
),
"Fairly High" = c(
100*round(phnorm(1, sigma = sd_w) - phnorm(0.5, sigma = sd_w),2)
),
"Fairly Extreme" = c(
100*round(1 - phnorm(1, sigma = sd_w), 2)
)) %>%
mutate(across(everything(), ~str_c(., "%"))) %>%
gt() %>%
cols_align(
align = "center",
columns = everything()
) %>%
# Add column tabs based on the column names defined above
tab_spanner(label = "Heterogeneity Range",
columns = c("Low",
"Reasonable",
"Fairly High",
"Fairly Extreme"))
| Heterogeneity Range | |||
|---|---|---|---|
| Low | Reasonable | Fairly High | Fairly Extreme |
| 16% | 52% | 27% | 5% |
# Fit the model
## Formula
mf =
formula(yi | se(sei) ~ 0 + oxygen + (1|study))
priors =
prior(normal(0, 0.82), class = "b") +
prior(normal(0, 0.5), class = "sd")
m1 =
brm(
data = d_logOR,
family = gaussian,
formula = mf,
prior = priors,
sample_prior = T,
control = list(adapt_delta = .97),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 100,
file = here("03_updated_analyses_who", "output", "fits", "m1.Rds"),
file_refit = "on_change"
)
# Detailed diagnostics
# launch_shinystan(m1)
Prior predictive check:
95% quantile intervals
m1 %>%
prior_draws() %>%
# Plot!
ggplot(aes(
# exp() to transform to Odds Ratio
x = exp(b))
) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(point_interval = median_qi,
.width = 0.95,
n = 1e4,
fill = "#6EB2E4") +
geom_vline(xintercept = 1, linetype = 2, size = 0.4, alpha = 0.7) +
labs(
x = "\nOdds Ratio",
y = "Density\n"
) +
scale_x_continuous(breaks = c(1, seq(from = 0, to = 5, 1))) +
coord_cartesian(xlim = c(0, 5.2)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 20, 20, 20))
diag_plot(model = m1,
pars_list = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV","sd_study__Intercept"),
ncol_trace = 4)
# Fit the model
## Formula
mf =
formula(yi | se(sei) ~ 0 + oxygen + (1|study))
priors =
prior(normal(0, 0.82), class = "b") +
prior(normal(0, 0.5), class = "sd")
m1_double =
brm(
data = d_logOR,
family = gaussian,
formula = mf,
prior = priors,
control = list(adapt_delta = .97),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 4000, # double
iter = 8000, # double
seed = 123,
file = here("03_updated_analyses_who", "output", "fits", "double_iterations",
"m1_double.Rds"),
file_refit = "on_change"
)
diag_plot(model = m1_double,
pars_list = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV", "sd_study__Intercept"),
ncol_trace = 4)
mcmc_hist(m1, pars = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV", "sd_study__Intercept"))
mcmc_acf(m1,
c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV", "sd_study__Intercept"),
lags = 10)
mcmc_dens(m1, pars = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV",
"sd_study__Intercept"),
transformations = list("b_oxygenlow" = "exp",
"b_oxygenNIV" = "exp",
"b_oxygenIMV" = "exp"))
Posterior predictive check
pp_check(m1, ndraws = 20)
Not applicable.
Here, we will fit models with different combinations of priors. For each parameter (\(\alpha\) and \(\tau\)), we will use three different types of prior:
Notably, we used weakly-informative priors for all of the parameters in our main model.
Here are the distributions used for \(\alpha\):
sd_w = 0.82
sd_v = 4
sd_i = 0.35
twosd_w = sd_w*1.96
twosd_v = sd_v*1.96
twosd_i = sd_i*1.96
# https://stackoverflow.com/questions/19950219/using-legend-with-stat-function-in-ggplot2
ggplot(data = data.frame(x = c(-5, 5)), aes(x)) + #Empty plot
# Normal(0, 0.82)
stat_function(fun = dnorm, n = 1000,
args = list(mean = 0, sd = sd_w), linetype=1, size = 1.2,
aes(colour = "0.82")) +
# Normal(0, 4)
stat_function(fun = dnorm, n = 1000,
args = list(mean = 0, sd = sd_v), linetype=1, size = 1.2,
aes(colour = "4.0")) +
# Normal(0, 0.35)
stat_function(fun = dnorm, n = 1000,
args = list(mean = 0, sd = sd_i), linetype=1, size = 1.2,
aes(colour = "0.35")) +
scale_colour_manual(latex2exp::TeX("Normal(0, $\\sigma)"),
values = pnw_palette("Starfish", 3, type = "discrete")) +
geom_vline(xintercept = c(-twosd_w, twosd_w),
linetype = 2,
color = pnw_palette("Starfish", 3, type = "discrete")[2]) +
geom_vline(xintercept = c(-twosd_i, twosd_i),
linetype = 2,
color = pnw_palette("Starfish", 3, type = "discrete")[1]) +
scale_y_continuous(breaks = seq(0, 1.5, 0.5),
limits = c(0, 1.5),
expand = c(0, 0)) + # remove gap between X and Y axis
scale_x_continuous(breaks = c(-3.92, -twosd_w, -twosd_i,
0, twosd_w, twosd_i, 3.92),
labels = function(x) round(as.numeric(x), 2),
expand = c(0, 0)) +
coord_cartesian(x = c(-5, 5)) +
labs(x = latex2exp::TeX("\n $\\alpha"),
y = "Density") +
theme_classic() +
theme(
plot.margin = margin(20,20,0,20),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 15),
legend.position = 'right',
legend.text = element_text(size=12),
legend.title = element_text(size=14),
legend.key= element_blank(),
panel.background = element_blank()
)
tibble(
"Prior" = c("Weakly-Informative", "Vague", "Informative"),
Mean = 0,
"SD" = c(sd_w, sd_v, sd_i)
) %>%
mutate("Mean " = exp(Mean),
"95% CI" =
str_c(" [",
round(exp(Mean - 1.96*`SD`), 1),
", ", round(exp(Mean + 1.96*`SD`), 1),
"]")) %>%
gt() %>%
cols_align(
align = "center",
columns = everything()
) %>%
cols_align(
align = "left",
columns = "Prior"
) %>%
# Add column tabs based on the column names defined above
tab_spanner(label = "Log Odds Ratio",
columns = c("Mean", "SD")) %>%
tab_spanner(label = "Odds Ratio",
columns = c("Mean ", "95% CI"))
| Prior | Log Odds Ratio | Odds Ratio | ||
|---|---|---|---|---|
| Mean | SD | Mean | 95% CI | |
| Weakly-Informative | 0 | 0.82 | 1 | [0.2, 5] |
| Vague | 0 | 4.00 | 1 | [0, 2540.2] |
| Informative | 0 | 0.35 | 1 | [0.5, 2] |
Here are the distributions used for \(\tau\):
sd_w = 0.5
sd_v = 4
informative = TurnerEtAlPrior("all-cause mortality",
"pharma", "placebo / control")
logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]
mean_i = -1.975 # logmean
sd_i = 0.67 # logsd
# https://stackoverflow.com/questions/19950219/using-legend-with-stat-function-in-ggplot2
ggplot(data = data.frame(x = c(0, 2)), aes(x)) + #Empty plot
# Half-Normal(0.5)
stat_function(fun = extraDistr::dhnorm, n = 1000,
args = list(sigma = sd_w), linetype=1, size = 1.2,
aes(colour = "Weakly informative")) +
# Half-Normal(4)
stat_function(fun = extraDistr::dhnorm, n = 1000,
args = list(sigma = sd_v), linetype=1, size = 1.2,
aes(colour = "Vague")) +
# Log-Normal(-1.975, 0.67)
stat_function(fun = dlnorm, n = 1000,
args = list(meanlog = mean_i,
sdlog = sd_i),
linetype=1, size = 1.2,
aes(colour = "Informative")) +
scale_colour_manual("Prior",
values = pnw_palette("Bay", 3, type = "continuous")) +
scale_y_continuous(breaks = seq(0, 6, 1.5),
limits = c(0, 6),
expand = c(0, 0)) + # remove gap between X and Y axis
scale_x_continuous(breaks = c(0, 0.1, 0.5, 1, 2),
labels = function(x) round(as.numeric(x), 1),
expand = c(0, 0)) +
coord_cartesian(x = c(0, 2)) +
labs(x = latex2exp::TeX("Between-study standard deviation ($\\tau)"),
y = "Density") +
theme_classic() +
theme(
plot.margin = margin(20,20,0,20),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 15),
legend.position = 'right',
legend.text = element_text(size=12),
legend.title = element_text(size=14),
legend.key= element_blank(),
panel.background = element_blank()
)
tibble(
"Prior" = c("Weakly-Informative", "Vague", "Informative"),
"Low" = c(
100*round(phnorm(0.1, sigma = sd_w) - phnorm(0, sigma = sd_w),2),
100*round(phnorm(0.1, sigma = sd_v) - phnorm(0, sigma = sd_v),2),
100*round(plnorm(0.1, meanlog = mean_i, sdlog = sd_i) -
plnorm(0, meanlog = mean_i, sdlog = sd_i),2)
),
"Reasonable" = c(
100*round(phnorm(0.5, sigma = sd_w) - phnorm(0.1, sigma = sd_w),2),
100*round(phnorm(0.5, sigma = sd_v) - phnorm(0.1, sigma = sd_v),2),
100*round(plnorm(0.5, meanlog = mean_i, sdlog = sd_i) -
plnorm(0.1, meanlog = mean_i, sdlog = sd_i),2)
),
"Fairly High" = c(
100*round(phnorm(1, sigma = sd_w) - phnorm(0.5, sigma = sd_w),2),
100*round(phnorm(1, sigma = sd_v) - phnorm(0.5, sigma = sd_v),2),
100*round(plnorm(1, meanlog = mean_i, sdlog = sd_i) -
plnorm(0.5, meanlog = mean_i, sdlog = sd_i),2)
),
"Fairly Extreme" = c(
100*round(1 - phnorm(1, sigma = sd_w), 2) ,
100*round(1 - phnorm(1, sigma = sd_v), 2),
100*round(1 - plnorm(1, meanlog = mean_i, sdlog = sd_i),2))
) %>%
mutate(across(2:last_col(), ~str_c(., "%"))) %>%
gt() %>%
cols_align(
align = "center",
columns = everything()
) %>%
cols_align(
align = "left",
columns = "Prior"
) %>%
# Add column tabs based on the column names defined above
tab_spanner(label = "Heterogeneity Range",
columns = c("Low",
"Reasonable",
"Fairly High",
"Fairly Extreme"))
| Prior | Heterogeneity Range | |||
|---|---|---|---|---|
| Low | Reasonable | Fairly High | Fairly Extreme | |
| Weakly-Informative | 16% | 52% | 27% | 5% |
| Vague | 2% | 8% | 10% | 80% |
| Informative | 31% | 66% | 3% | 0% |
As mentioned above, we used all weakly-informative priors for our main model (already fitted). Now, we will fit two more models:
## Priors
# alpha
alpha_vague = prior(normal(0, 4), class = "b")
alpha_informative = prior(normal(0, 0.35), class = "b")
# Between-study standard deviation (tau_study)
# the package bayesmeta has a handy function to get the informative prior
informative = TurnerEtAlPrior("all-cause mortality", "pharma", "placebo / control")
logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]
tau_vague = prior(normal(0, 4), class = "sd", group = "study")
tau_informative = prior(lognormal(-1.975, # logmean
0.67), # logsd
class = "sd", group = "study")
# Put them all together in every possible combination, 1 line per model
priors = list(
# All vague
c(alpha_vague, tau_vague),
# All informative
c(alpha_informative, tau_informative)
)
### Create list to store fits
sensitivity_models = list()
#
### Store the main model
#
sensitivity_models[[1]] = m1
##
labs = data.frame(x = c("Weakly informative", "Vague", "Informative"))
# #
### Supressed because we saved the models already
#
# # Run the loop, to fit other models with different priors
#
# for (i in 1:(nrow(labs) - 1)) {
#
# # Show what model is running
# print(i)
#
# # Skip first index because main model (m1) is already stored there
# k = i + 1
#
# # Re-fit the main model, but now changing the prior
# new_fit = update(m1,
# prior = priors[[i]],
#
# cores = parallel::detectCores(),
# chains = 4,
# control = list(adapt_delta = .99),
# backend = "cmdstanr",
# seed = 123)
#
# # Save models
# sensitivity_models[[k]] = new_fit
# }
#
# names(sensitivity_models) = labs$x
#
# saveRDS(sensitivity_models,
# here("03_updated_analyses_who", "output", "fits", "sensitivity.Rds"))
sensitivity_models =
readRDS(here("03_updated_analyses_who", "output", "fits", "sensitivity.Rds"))
P.S.: The diagnostics plots of these two models will be shown in the end of this document.
Let’s visualize the prior and posterior distributions (constructed by the samples stored in each model object).
weakly_prior =
sensitivity_models[["Weakly informative"]] %>%
prior_draws() %>%
select(b) %>%
rename("Weakly informative" = b)
vague_prior =
sensitivity_models[["Vague"]] %>%
prior_draws() %>%
select(b) %>%
rename("Vague" = b)
informative_prior =
sensitivity_models[["Informative"]] %>%
prior_draws() %>%
select(b) %>%
rename("Informative" = b)
priors =
bind_cols(weakly_prior, vague_prior, informative_prior) %>%
pivot_longer(1:3) %>%
ggplot(aes(value)) +
stat_halfeye(point_interval = median_hdi, .width = .95,
fill = "#347178", slab_color = "grey92",
breaks = 10, slab_size = .2, outline_bars = T) +
scale_x_continuous(breaks = seq(-3, 3, 1)) +
scale_y_continuous(breaks = seq(0,1, 0.5)) +
coord_cartesian(x = c(-3, 3)) +
labs(x = latex2exp::TeX("$\\alpha"),
y = "Density\n") +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 20, 20, 20)
) +
facet_wrap(~name, ncol = 1)
weakly_posterior =
sensitivity_models[["Weakly informative"]] %>%
tidy_draws() %>%
select(b_oxygenlow, b_oxygenNIV, b_oxygenIMV) %>%
rename("Weakly informative [k = 1]" = b_oxygenlow,
"Weakly informative [k = 2]" = b_oxygenNIV,
"Weakly informative [k = 3]" = b_oxygenIMV)
vague_posterior =
sensitivity_models[["Vague"]] %>%
tidy_draws() %>%
select(b_oxygenlow, b_oxygenNIV, b_oxygenIMV) %>%
rename("Vague [k = 1]" = b_oxygenlow,
"Vague [k = 2]" = b_oxygenNIV,
"Vague [k = 3]" = b_oxygenIMV)
informative_posterior =
sensitivity_models[["Informative"]] %>%
tidy_draws() %>%
select(b_oxygenlow, b_oxygenNIV, b_oxygenIMV) %>%
rename("Informative [k = 1]" = b_oxygenlow,
"Informative [k = 2]" = b_oxygenNIV,
"Informative [k = 3]" = b_oxygenIMV)
posteriors =
bind_cols(weakly_posterior, vague_posterior, informative_posterior) %>%
pivot_longer(1:last_col()) %>%
ggplot(aes(value)) +
stat_halfeye(point_interval = median_hdi, .width = .95,
fill = "#605A91", slab_color = "grey92",
breaks = 10, slab_size = .2, outline_bars = T) +
scale_x_continuous(breaks = seq(-1, 1, 0.5)) +
scale_y_continuous(breaks = seq(0,1, 0.5)) +
coord_cartesian(x = c(-1, 1)) +
labs(x = latex2exp::TeX("$\\alpha"),
y = "Density\n") +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 20, 20, 20)
) +
facet_wrap(~name, ncol = 3)
priors + posteriors +
plot_annotation(tag_levels = "A")
Panel A: Prior distributions; Panel B: Posterior distributions; k = 1 represents the simple oxygen only subgroup, k = 2 noninvasive ventilation, k = 3 invasive mechanical ventilation
weakly_prior =
sensitivity_models[["Weakly informative"]] %>%
prior_draws() %>%
select(sd_study) %>%
rename("Weakly informative" = sd_study)
vague_prior =
sensitivity_models[["Vague"]] %>%
prior_draws() %>%
select(sd_study) %>%
rename("Vague" = sd_study)
vague_prior_hdi =
vague_prior %>%
median_hdi()
# 2.62 [0.000237, 7.77]
informative_prior =
sensitivity_models[["Informative"]] %>%
prior_draws() %>%
select(sd_study) %>%
rename("Informative" = sd_study)
priors =
bind_cols(weakly_prior, vague_prior, informative_prior) %>%
pivot_longer(1:3) %>%
mutate(name = fct_rev(name)) %>%
ggplot(aes(value)) +
stat_halfeye(point_interval = median_hdi, .width = .95,
fill = "#9B5446", slab_color = "grey92",
breaks = 10, slab_size = .2, outline_bars = T) +
scale_x_continuous(breaks = seq(0, 1, 0.25)) +
scale_y_continuous(breaks = seq(0,1, 0.5)) +
coord_cartesian(x = c(0, 1)) +
labs(x = latex2exp::TeX("Between-study standard deviation ($\\tau)"),
y = "Density\n") +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 20, 20, 20)
) +
facet_wrap(~name, ncol = 1)
weakly_posterior =
sensitivity_models[["Weakly informative"]] %>%
tidy_draws() %>%
select(sd_study__Intercept) %>%
rename("Weakly informative" = sd_study__Intercept)
vague_posterior =
sensitivity_models[["Vague"]] %>%
tidy_draws() %>%
select(sd_study__Intercept) %>%
rename("Vague" = sd_study__Intercept)
informative_posterior =
sensitivity_models[["Informative"]] %>%
tidy_draws() %>%
select(sd_study__Intercept) %>%
rename("Informative" = sd_study__Intercept)
posteriors =
bind_cols(weakly_posterior, vague_posterior, informative_posterior) %>%
pivot_longer(1:3) %>%
mutate(name = fct_rev(name)) %>%
ggplot(aes(value)) +
stat_halfeye(point_interval = median_hdi, .width = .95,
fill = "#9B5446", slab_color = "grey92",
breaks = 10, slab_size = .2, outline_bars = T) +
scale_x_continuous(breaks = seq(0, 1, 0.25)) +
scale_y_continuous(breaks = seq(0,1, 0.5)) +
coord_cartesian(x = c(0, 1)) +
labs(x = latex2exp::TeX("Between-study standard deviation ($\\tau)"),
y = "Density\n") +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'none',
plot.margin = margin(20, 20, 20, 20)
) +
facet_wrap(~name, ncol = 1)
priors + posteriors +
plot_annotation(tag_levels = "A")
Panel A: Prior distributions; Panel B: Posterior distributions
We did not fit models “adjusting the entire prior distribution (i.e., using a completely different prior distribution than before) or adjusting hyperparameters upward and downward and reestimating the model with these varied priors.”
On the models shown in “Point 8: Is There a Notable Effect of the Prior When Compared With Noninformative Priors?”
# Run loop
for (i in 1:nrow(labs)) {
# Run custom function diag_plot, 1 per model
p = diag_plot(model = sensitivity_models[[i]],
pars_list = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV", "sd_study__Intercept"),
ncol_trace = 4)
# Display plot
print(p)
# Add legend to show respective priors
lab = paste0(labs %>% slice(i) %>% pull() )
grid::grid.text(lab, 0.25, .03, gp=grid::gpar(cex=1.3))
}
for (i in 1:nrow(labs)) {
p = mcmc_hist(sensitivity_models[[i]],
pars = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV", "sd_study__Intercept"))
# Display plot
print(p)
lab = paste0(labs %>% slice(i) %>% pull() )
grid::grid.text(lab, 0.82, .2, gp=grid::gpar(cex=1.3))
}
Posterior predictive check
# Run lopp
for (i in 1:nrow(labs)) {
p = pp_check(sensitivity_models[[i]], ndraws = 20)
print(p)
lab = paste0(labs %>% slice(i) %>% pull() )
grid::grid.text(lab, .7, .9, gp=grid::gpar(cex=1.2))
}
for (i in 1:nrow(labs)) {
p = mcmc_acf(sensitivity_models[[i]],
pars = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV",
"sd_study__Intercept"),
lags = 10)
# Display plot
print(p)
lab = paste0(labs %>% slice(i) %>% pull() )
grid::grid.text(lab, 0.25, .03, gp=grid::gpar(cex=1.3))
}
for (i in 1:nrow(labs)) {
p = mcmc_dens(m1, pars = c("b_oxygenlow", "b_oxygenNIV", "b_oxygenIMV",
"sd_study__Intercept"),
transformations = list("b_oxygenlow" = "exp",
"b_oxygenNIV" = "exp",
"b_oxygenIMV" = "exp"))
# Display plot
print(p)
lab = paste0(labs %>% slice(i) %>% pull() )
grid::grid.text(lab, 0.82, .2, gp=grid::gpar(cex=1.3))
}
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices datasets utils methods
## [8] base
##
## other attached packages:
## [1] bayesplot_1.8.1 cowplot_1.1.1 Hmisc_4.5-0
## [4] Formula_1.2-4 survival_3.2-12 lattice_0.20-44
## [7] latex2exp_0.5.0 extraDistr_1.9.1 tidybayes_3.0.1
## [10] patchwork_1.1.1 gt_0.3.1 PNWColors_0.1.0
## [13] bayesmeta_2.7 numDeriv_2016.8-1.1 metafor_3.0-2
## [16] Matrix_1.3-4 forestplot_2.0.1 checkmate_2.0.0
## [19] magrittr_2.0.1 rio_0.5.27 here_1.0.1
## [22] brms_2.16.1 Rcpp_1.0.7 forcats_0.5.1
## [25] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
## [28] readr_2.0.1 tidyr_1.1.3 tibble_3.1.4
## [31] ggplot2_3.3.5 tidyverse_1.3.1 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.1 lme4_1.1-27.1
## [4] htmlwidgets_1.5.4 munsell_0.5.0 codetools_0.2-18
## [7] DT_0.19 miniUI_0.1.1.1 withr_2.4.2
## [10] Brobdingnag_1.2-6 colorspace_2.0-2 highr_0.9
## [13] knitr_1.34 rstudioapi_0.13 stats4_4.1.2
## [16] labeling_0.4.2 rstan_2.21.2 farver_2.1.0
## [19] bridgesampling_1.1-2 rprojroot_2.0.2 coda_0.19-4
## [22] vctrs_0.3.8 generics_0.1.0 xfun_0.25
## [25] R6_2.5.1 markdown_1.1 HDInterval_0.2.2
## [28] gamm4_0.2-6 projpred_2.0.2 assertthat_0.2.1
## [31] promises_1.2.0.1 scales_1.1.1 nnet_7.3-16
## [34] gtable_0.3.0 processx_3.5.2 rlang_0.4.11
## [37] splines_4.1.2 broom_0.7.9 inline_0.3.19
## [40] yaml_2.2.1 reshape2_1.4.4 abind_1.4-5
## [43] modelr_0.1.8 threejs_0.3.3 crosstalk_1.1.1
## [46] backports_1.2.1 httpuv_1.6.3 rsconnect_0.8.24
## [49] tensorA_0.36.2 tools_4.1.2 ellipsis_0.3.2
## [52] jquerylib_0.1.4 posterior_1.1.0 RColorBrewer_1.1-2
## [55] ggridges_0.5.3 plyr_1.8.6 base64enc_0.1-3
## [58] ps_1.6.0 prettyunits_1.1.1 rpart_4.1-15
## [61] zoo_1.8-9 haven_2.4.3 cluster_2.1.2
## [64] fs_1.5.0 data.table_1.14.0 ggdist_3.0.0
## [67] openxlsx_4.2.4 colourpicker_1.1.0 reprex_2.0.1
## [70] mvtnorm_1.1-2 matrixStats_0.60.1 hms_1.1.0
## [73] shinyjs_2.0.0 mime_0.11 evaluate_0.14
## [76] arrayhelpers_1.1-0 xtable_1.8-4 shinystan_2.5.0
## [79] jpeg_0.1-9 readxl_1.3.1 gridExtra_2.3
## [82] rstantools_2.1.1 compiler_4.1.2 V8_3.4.2
## [85] crayon_1.4.1 minqa_1.2.4 StanHeaders_2.21.0-7
## [88] htmltools_0.5.2 mgcv_1.8-36 later_1.3.0
## [91] tzdb_0.1.2 RcppParallel_5.1.4 lubridate_1.7.10
## [94] DBI_1.1.1 dbplyr_2.1.1 MASS_7.3-54
## [97] boot_1.3-28 cli_3.0.1 parallel_4.1.2
## [100] igraph_1.2.6 pkgconfig_2.0.3 foreign_0.8-81
## [103] xml2_1.3.2 svUnit_1.0.6 dygraphs_1.1.1.6
## [106] bslib_0.3.0 rvest_1.0.1 distributional_0.2.2
## [109] callr_3.7.0 digest_0.6.27 rmarkdown_2.10
## [112] cellranger_1.1.0 htmlTable_2.2.1 curl_4.3.2
## [115] shiny_1.6.0 gtools_3.9.2 nloptr_1.2.2.2
## [118] lifecycle_1.0.0 nlme_3.1-152 jsonlite_1.7.2
## [121] cmdstanr_0.4.0.9001 fansi_0.5.0 pillar_1.6.2
## [124] loo_2.4.1 fastmap_1.1.0 httr_1.4.2
## [127] pkgbuild_1.2.0 glue_1.4.2 xts_0.12.1
## [130] zip_2.2.0 png_0.1-7 shinythemes_1.2.0
## [133] stringi_1.7.4 sass_0.4.0 latticeExtra_0.6-29
## [136] renv_0.14.0 mathjaxr_1.4-0